GIS 与密铺与铺瓷砖的紧密关系
| 前言：铺瓷砖是一门古老的技术，与我们息息相关，其中蕴含的数学道理更与 GIS 紧密结合... |
什么是密铺？
密铺，什么是密铺，可以简单理解成铺瓷砖。
就是在一个平面上不留缝隙，把平面全部铺上，虽然咋一听上去不是什么大事，但是密铺是一门非常古老的研究，从古希腊就有记载（可能当时也是研究铺地板怎么好看吧），同时密铺在数学界有着严谨的分类和定义。
我们这里说的密铺包括下文的结论均指平面多边形单密铺。
平面多边形单密铺指使用一种凸多边形图形实现密铺，直接说一些基本结论吧：
任何三角形都可以密铺整个平面；
任何凸四边形（包括正方形，矩形）都可以密铺整个平面；
正五边形不能密铺；
只有正六边形能密铺平面；
任何凸 n 边形（n>=7)都没有合适的单密铺方式；
...
所以只需要知道有多少种五边形能密铺，就宣告人类发现了所有的多边形单密铺，这个历程跨越百年；直到2015年，华盛顿大学的教授 Casey Mann 及其妻子 Jennifer McLoud，再加上学生 David Von Derau 发现第15种可以平面密铺的五边形。
尽管如此，人们还是不清楚能够密铺的五边形到底还有没有，毕竟这次发现距离上次已经过了30年，没有人能保证下一个30年不会发现第16种。
直达2017年，法国数学家 Michaël Rao 通过计算机穷举最后得出了结论并证明，那就是能够平面密铺的五边形只有这15种。
自此，人类已经完全发现了能够实现平面多边形单密铺的所有图形。
密铺与 GIS
与GIS 的关系
那么和 GIS 有什么关系呢？
密铺和 GIS 的关系太多了，我们平时就有很多使用到密铺的地方，比如 创建渔网、格网 等，简单来说就是方形格网，就像下面这样，常见的还有六边形网格。
这些网格作用很大，可用于很多方面，比如栅格数据按指定形状重采样、区域划分等等。
Python 制作密铺形状
ArcGIS 可以创建矩形的密铺网格，10.3之后的版本可以创建六边形密铺网格。
QGIS 3.16 版本除了可以创建矩形、六边形外，还可以创建菱形密铺。
但是呢，我还没有发现五边形密铺，所以我使用 Python 写了两种密铺五边形的创建脚本，结合 ArcPy 可以创建矢量文件，效果如下：
如何创建这种形状的几何呢？我们以第一个形状为例：
首先分析单一图形，上面两条长边是下面短边长的两倍，最上面的角是60°，剩下的都是120°。知道了这些基础信息，就能在坐标系中准确定位了，只有在坐标系中确定了点的位置，才能继而构建线或者面几何。
坐标系使用的是笛卡尔坐标系（就是最常见的数学坐标系），向上为 Y 轴正方向，向右为 X 轴正方向。在这个坐标系下，现在可以轻松的计算出A、B、C、D、E点（示意图1）的坐标
from __future__ import division
from math import sin, radians, pow, sqrt
# 短边
length = 300
# 角度
angle01 = 60
# 垂直距离
leng = length * sin(radians(angle01)) # 300*sin60°
pta = (oX + length * 2, oY)
ptb = (pta[0] + length / 2, oY + leng)
ptc = (pta[0], oY + leng * 2)
ptd = (oX + length, ptc[1])
pte = (oX + length / 2, oY + leng * 3)
其中 oX、oY 表示初始点的 x、y 坐标值；
pta 表示 A点（示意图1）坐标，后面依次对应。
当然不可能每个五边形的点都算一次，找到一个有规律的整体，然后平移再平移就行了。
示意图1中的1、2、3、4、5、6块五边形组成了一个整体，我们可以轻松的获取到这六个五边形各个点的坐标，然后直接平移整体即可完美覆盖所有区域。
所以现在搞清楚平移距离和方向以及平移次数就完成了。
o 点到 M 点（示意图1绿色线），即整体向正 Y 方向平移，o 点到 N 点（示意图1绿色线），即整体向正 X 方向平移；这种平移可以使用 numpy 进行坐标加减，效率很高，所以这两个脚本的运行效率几乎可以媲美 ArcGIS 的原生工具。
好的，基本思路就是这样，下面上完整代码，感兴趣的读者可以仔细看看，有空我会更新之前的样式箱，加入这几种创建五边形密铺的小工具，到时候会发布更新文章：
# -------------------------------------------
# Name: tile
# Author: Hygnic
# Created on: 2021/7/22 11:11
# Version:
# Reference:
"""
Description:
Usage:
"""
# -------------------------------------------
from __future__ import absolute_import
from __future__ import division
import os
import arcpy
import numpy as np
from math import sin, radians, pow, sqrt
"""--------------------------------------"""
"""--------基本方法------"""
def check_extent(input):
"""
获取输入图层本身的尺寸
:param input: 输入图层地址
:return: 原点，宽，高
"""
lyr = arcpy.mapping.Layer(input)
lyr_extent = lyr.getExtent()
_origin = (lyr_extent.XMin, lyr_extent.YMin)
if lyr.isRasterLayer:
# arcpy.env.extent = lyr_extent
# arcpy.env.mask = lyr
pass
return _origin, lyr_extent.width, \
lyr_extent.height, lyr_extent.spatialReference
def tile_creator(array_obj, featurecalss):
"""
将 array_obj 中的点去除然后生成面
:param array_obj:
:param featurecalss: 矢量文件
:return:
"""
_rows = arcpy.da.InsertCursor(featurecalss, "SHAPE@")
for _ii in array_obj:
_rows.insertRow([_ii])
del _rows
"""--------基本方法------"""
"""--------------------------------------"""
"""--------------------------------------"""
"""--------基本属性------"""
ws = os.path.abspath(os.getcwd())
arcpy.env.workspace = ws
arcpy.env.overwriteOutput = True
# 坐标原点
# lyr_o, lyr_w, lyr_h, sr=check_extent("data/grid.shp")
lyr_o, lyr_w, lyr_h, sr=check_extent("tiff.tif")
print "featureclass width:{}".format(lyr_w)
print "featureclass height:{}".format(lyr_h)
origin = lyr_o
oX = origin[0]
oY = origin[1]
print "origin X:{}".format(oX)
print "origin Y:{}".format(oY)
# 五边形短边长度
length = 300
# 角度
angle01 = 60
# 垂直距离
leng = length * sin(radians(angle01)) # 300*sin60°
"""--------基本属性------"""
"""--------------------------------------"""
"""--------------------------------------"""
"""--------基本坐标------"""
#一组五边形在坐标四个象限内的坐标
pta = (oX + length * 2, oY)
ptb = (pta[0] + length / 2, oY + leng)
ptc = (pta[0], oY + leng * 2)
ptd = (oX + length, ptc[1])
pte = (oX + length / 2, oY + leng * 3)
quadrant1 = [origin, pta, ptb, ptc, ptd, pte]
# 二象限
pta2 = (oX - length * 2, pta[1])
ptb2 = (pta2[0] - length / 2, ptb[1])
ptc2 = (pta2[0], ptc[1])
ptd2 = (oX - length, ptd[1])
pte2 = (oX - length / 2, pte[1])
quadrant2 = [origin, pta2, ptb2, ptc2, ptd2, pte2]
# 三象限
ptb3 = (ptb2[0], oY - leng)
ptc3 = (ptc2[0], oY - leng * 2)
ptd3 = (ptd2[0], ptc3[1])
pte3 = (pte2[0], oY - leng * 3)
quadrant3 = [origin, pta2, ptb3, ptc3, ptd3, pte3]
# 四象限
ptb4 = (pta[0] + length / 2, oY - leng)
ptc4 = (pta[0], oY - leng * 2)
ptd4 = (oX + length, ptc4[1])
pte4 = (oX + length / 2, oY - leng * 3)
quadrant4 = [origin, pta, ptb4, ptc4, ptd4, pte4]
pts = [origin, pta, ptb, ptc, ptd, pte,
origin, pta2, ptb2, ptc2, ptd2, pte2,
origin, pta2, ptb3, ptc3, ptd3, pte3,
origin, pta, ptb4, ptc4, ptd4, pte4]
"""--------基本坐标------"""
"""--------------------------------------"""
"""--------------------------------------"""
"""--------构建要素------"""
# 重要参数和属性
cfm = arcpy.CreateFeatureclass_management
shpfile = cfm(ws, "PentagonTile", "polygon", spatial_reference=sr)
# 左上方向的偏移距离
offset_x = -length * 3/2
offset_y = 5 * leng
# 右下方向的偏移距离
offset_x2 = length*4.5
offset_y2 = -leng
# y轴方向循环次数
loop_y = int(lyr_h/(6*leng))
array_pt = [] # 用于存放一整列的五边形
for i in xrange(int(loop_y*1.6)):
# 向上偏移距离
# new_pts = [(_[0] - length * 3 / 2 * i, _[1] + 5 * leng * i) for _ in pts]
new_pts = [(_[0]+offset_x*i, _[1]+offset_y*i) for _ in pts]
A = new_pts[:5]
B = new_pts[4:6] + [new_pts[11], new_pts[10], new_pts[0]]
C = new_pts[6:11]
D = new_pts[12:17]
E = new_pts[16:18] + [new_pts[-1], new_pts[-2], new_pts[0]]
F = new_pts[-6:-1]
array_pt.append(A)
array_pt.append(B)
array_pt.append(C)
array_pt.append(D)
array_pt.append(E)
array_pt.append(F)
np_array = np.array(array_pt)
print "Len:{}".format(len(np_array)) # 396
print "Size:{}".format(np_array.size)
print "Shape:{}".format(np_array.shape)
print "Ndim:{}".format(np_array.ndim)
loop_x = int(lyr_w/(4*length))
for _ in xrange(int(loop_x*1.4)):
tile_creator(np_array, shpfile)
np_array = np_array+(offset_x2, offset_y2)
"""--------构建要素------"""
"""--------------------------------------"""
结束语
这篇文章的代码有三点价值：
使用 ArcPy 创建几何；
使用 numpy 做简单的矩阵运算；
实现了五边形的密铺，提供了更多选择，不仅仅局限于矩形等。
老环节，随机推荐以往文章：
《ArcPy教程指南2021》——从0开始，从实际应用出发，带大家熟悉掌握 ArcPy。（持续更新中！）
参考：
[1]Pentagon Tiling Proof Solves Century-Old Math Problem.https://www.quantamagazine.org/pentagon-tiling-proof-solves-century-old-math-problem-20170711/
[2]Tiling-Github.https://github.com/hygnic/Tiling
[3]如何看待卡西·曼夫妇发现的可无缝密铺平面的五边形？https://www.zhihu.com/question/34916541/answer/60644653
[4]Cairo Tilings with minimal math.https://medium.com/@femion/cairo-tilings-with-minimal-math-bee322e716e1
[5]Pentagonal tiling wiki.https://en.wikipedia.org/wiki/Pentagonal_tiling